Non-linear models: interactive

PLS 206 — Applied Multivariate Modeling

Three explorers for Week 4. They run in your browser with no R and no download, so they work on a phone too. Tick “Show R code” in any of them to get code that reproduces the current view in RStudio. The first two use the growth-chamber seedlings from Monday (seedling_temp.csv, standard variety); the third uses the soil collars (soil_season.csv).

Cheat sheet: reading summary() of a GAM

The numbers are from gam(mass_mg ~ s(temp), data = std, method = "REML") on the 150 standard-variety pots.

line of summary(m) here what it means
Family: gaussian, Link function: identity the GLM part: the distribution of \(y\) and the link, exactly as in Week 3
(Intercept) 42.68 the mean of \(y\): every smooth is centred to average zero, so the intercept is the overall mean shoot mass
edf 6.87 effective degrees of freedom: how many parameters the curve is worth after the penalty; 1 = a straight line
Ref.df 7.98 reference degrees of freedom used for the test; not something to interpret
F, p-value 192.4, < 2e-16 an approximate test of “is the smooth flat?” (no effect of temperature at all)
R-sq.(adj) 0.912 adjusted R², as in lm()
Deviance explained 91.6% the share of the null deviance the model explains (for a Normal GAM, close to R²)
-REML 450.4 the score REML minimised to choose λ; only useful for comparing fits of the same data
Scale est. 21.26 the residual variance \(\hat\sigma^2\) (here in mg²; its square root, 4.6 mg, is the residual SD)
n 150 number of observations

Is k big enough? k is a ceiling on wiggliness, not the wiggliness: the penalty decides how much of it is used.

  1. Start generous (the default k = 10 is usually fine for one predictor) and run gam.check(m) or k.check(m).
  2. Worry only if the k-index is below 1 with a small p-value and edf is close to k' (that is, k − 1).
  3. Then double k and refit. If the curve barely changes, you were fine. Here k = 10 passes (k-index 0.95, p well above 0.05; the p-value comes from a random shuffle, so it moves a little from run to run); k = 4 fails (k-index 0.79, p = 0.0025, edf 2.99 of a possible 3).

Bumps, weights and the penalty

A spline is k bumps, each given a weight, added up. The thin coloured curves are the bumps after weighting; the thick blue curve is their sum, which is the fitted curve. The penalty λ is the price of a wiggle: move it and watch the weights shrink until the curve is a straight line. Switch “Who sets the penalty λ?” to REML to let the data choose, as gam(..., method = "REML") does, and then change k: once k is generous, the curve barely moves.

The edf readout is the trace of the hat matrix minus one (the intercept is counted separately, as in mgcv’s summary). The small plot underneath shows edf for every λ, with the current λ in red. The bumps here are P-splines (mgcv’s bs = "ps"), so the numbers differ a little from the default s(temp) in the slides, which uses a different set of bumps.

Data

Who sets the penalty λ?

Global polynomials, local splines

Both models are lm() with extra columns. On the left the columns are powers of temperature, poly(temp, degree, raw = TRUE); on the right they are spline bumps, bs(temp, df). The bottom panels show those columns: every power of temperature is non-zero across the whole range, while each spline bump only lives in one neighbourhood. That difference shows up in two places. Past the edge of the data (the shaded band is where the pots are) a high-degree polynomial swings off to absurd values. And if you pick 30 pots and nudge the coldest one, the polynomial changes at the hot end too, while the spline only changes near the pot you moved. (Past the data, bs() carries on along its last cubic piece, so it drifts too, just far less; a natural spline, ns(), carries on as a straight line instead.)

Pots

The faint dashed line is the fit before the nudge.

Cycles: a curve that meets itself

Day of year goes round in a circle: 31 December sits next to 1 January. An ordinary spline does not know that, so its two ends are free to disagree. A cyclic spline uses bumps that wrap from December into January, so the fitted curve always joins up. “Wrap the year” repeats the data on either side of the year (faint points) so you can see the seam. λ is chosen by REML for both fits. These are cyclic P-splines (bs = "cp"); the lecture uses bs = "cc", which works the same way with a slightly different set of bumps. The model here ignores the eight collars; the lecture adds them with s(collar, bs = "re").

λ is chosen by REML for both.

Show the bumps of